Tutorial Part III — SDM with MaxEnt

Course 3 SDM with MaxEnt

Preparations

Load libraries

### The packages (repetition) ###
   library(here)
   library(sp)
   library(sf)
   library(tmap)
   library(dismo)
   library(raster)
   library(GISTools)
   library(rgdal)
   library(maptools)
   library(viridis)
   library(rgeos)
   library(rgl)
 # library(rasterVis) 
   library(unmarked)
   library(spatstat)

## check that the package terra is detached!
# detach(package:terra)

Set the workspace

### The workspace  (repetition) ###
getwd() # you can also use the package 'here()'
## [1] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/Course3_SDM/R"
# my data are outside my course-folder, therefore I have to do it the old-fashioned way.
root_wd <- "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub"
# relative to work-wd
maps_wd <- paste(root_wd,"/","_DataBorneo/BaseMaps",sep='')
recs_wd <- paste(root_wd,"/","_DataBorneo/BaseRecords",sep='')

# the output folder should have been created by you during Tutorial 2 'R goes spatial'.
# It should contain the hillshade.asc
output_wd <- paste(root_wd,"/","output",sep='')

setwd(output_wd) ## set to the OUTPUT folder!
getwd()
## [1] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/output"

Load spatial data and define the CRS

Load rasters

### Data import - Load spatial data (repetition) ###
## note, here I use the 'dialect' of the old 'raster'-package.

ras1 <- raster(paste(maps_wd,'/bio_asc_01.asc',sep=''))
# assign the projection (crs - coordinate reference system)
ras1@crs <- CRS("+proj=longlat +datum=WGS84")

ras24 <- raster(paste(maps_wd,'/bio_asc_24.asc',sep='')) #DEM
ras24@crs <- ras1@crs # same as: CRS("+proj=longlat +datum=WGS84")

ras42 <- raster(x = paste(maps_wd,'/bio_asc_42.asc',sep='')) 
ras42@crs <- CRS(projargs = "+init=epsg:4326") 

hillsh <- raster(x = paste(output_wd,'/hillshade.asc',sep='')) 
hillsh@crs <- CRS(projargs = "+init=epsg:4326")

Stack environmental variables

Loading many rasters can be done in one go.

### Raster stack of environmental predictors (repetition) ###

# the list.files command is very helpful to check what is in the folders
# use 'pattern' for searching for special file types, here .asc-files:
files <- list.files(path= maps_wd, pattern='.asc$',
                    full.names=TRUE )
files # these are just names! To load them as spatial objects, use raster() or stack()
##  [1] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/_DataBorneo/BaseMaps/bio_asc_01.asc"
##  [2] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/_DataBorneo/BaseMaps/bio_asc_02.asc"
##  [3] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/_DataBorneo/BaseMaps/bio_asc_03.asc"
##  [4] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/_DataBorneo/BaseMaps/bio_asc_04.asc"
##  [5] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/_DataBorneo/BaseMaps/bio_asc_05.asc"
##  [6] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/_DataBorneo/BaseMaps/bio_asc_06.asc"
##  [7] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/_DataBorneo/BaseMaps/bio_asc_07.asc"
##  [8] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/_DataBorneo/BaseMaps/bio_asc_08.asc"
##  [9] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/_DataBorneo/BaseMaps/bio_asc_09.asc"
## [10] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/_DataBorneo/BaseMaps/bio_asc_10.asc"
## [11] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/_DataBorneo/BaseMaps/bio_asc_11.asc"
## [12] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/_DataBorneo/BaseMaps/bio_asc_12.asc"
## [13] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/_DataBorneo/BaseMaps/bio_asc_13.asc"
## [14] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/_DataBorneo/BaseMaps/bio_asc_14.asc"
## [15] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/_DataBorneo/BaseMaps/bio_asc_15.asc"
## [16] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/_DataBorneo/BaseMaps/bio_asc_16.asc"
## [17] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/_DataBorneo/BaseMaps/bio_asc_17.asc"
## [18] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/_DataBorneo/BaseMaps/bio_asc_18.asc"
## [19] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/_DataBorneo/BaseMaps/bio_asc_19.asc"
## [20] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/_DataBorneo/BaseMaps/bio_asc_21.asc"
## [21] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/_DataBorneo/BaseMaps/bio_asc_22.asc"
## [22] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/_DataBorneo/BaseMaps/bio_asc_24.asc"
## [23] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/_DataBorneo/BaseMaps/bio_asc_27.asc"
## [24] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/_DataBorneo/BaseMaps/bio_asc_42.asc"
## [25] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/_DataBorneo/BaseMaps/hillshade.asc"
predictors <- stack(files[c(9,12,22,24)]) # note that I only load 4 rasters

# set CRS:
predictors@crs <- CRS(projargs = "+init=epsg:4326")

Plot rasters

plot(predictors, col= viridis(100)) # might take some time depending on your computer

# what are the differences? Have a look also at the data description file in the data folder.

Load shapefiles

### Read in some Shapefiles (repetition) ###

# package sp (old, but easier for plotting)
Borneo_shp <- readOGR(dsn=maps_wd, layer="borneo_admin",
              stringsAsFactors=FALSE)[,c(1:3,5,7,17,18)] #only load some of the columns of df
## OGR data source with driver: ESRI Shapefile 
## Source: "E:\PopDynIZW Dropbox\Steph Kramer\_GitHub\_DataBorneo\BaseMaps", layer: "borneo_admin"
## with 158 features
## It has 18 fields
PA_shp <-  readOGR(dsn=maps_wd, layer="Bor_PA",
                   stringsAsFactors=FALSE)[,c(1:4)]
## OGR data source with driver: ESRI Shapefile 
## Source: "E:\PopDynIZW Dropbox\Steph Kramer\_GitHub\_DataBorneo\BaseMaps", layer: "Bor_PA"
## with 255 features
## It has 12 fields
River_shp <- readOGR(dsn=maps_wd, layer="sn_100000",
                     stringsAsFactors=FALSE)
## OGR data source with driver: ESRI Shapefile 
## Source: "E:\PopDynIZW Dropbox\Steph Kramer\_GitHub\_DataBorneo\BaseMaps", layer: "sn_100000"
## with 396 features
## It has 4 fields
# package sf (new and maintained)
# Borneo outline polygon
Borneo_shp_sf <- st_read(dsn = maps_wd, 
                      layer = "borneo_admin",
                      stringsAsFactors = FALSE)[,c(1:3,5,7,17,18)]
## Reading layer `borneo_admin' from data source `E:\PopDynIZW Dropbox\Steph Kramer\_GitHub\_DataBorneo\BaseMaps' using driver `ESRI Shapefile'
## Simple feature collection with 158 features and 18 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 108.5986 ymin: -4.943054 xmax: 119.2697 ymax: 7.380556
## geographic CRS: WGS 84
# Protected areas (PA) polygon
PA_shp_sf <-  st_read(dsn = maps_wd, 
                   layer = "Bor_PA",
                   stringsAsFactors = FALSE)[, c(1:4)]
## Reading layer `Bor_PA' from data source `E:\PopDynIZW Dropbox\Steph Kramer\_GitHub\_DataBorneo\BaseMaps' using driver `ESRI Shapefile'
## Simple feature collection with 255 features and 12 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 108.5969 ymin: -4.18356 xmax: 119.6176 ymax: 9.911229
## geographic CRS: WGS 84
# Rivers lines
River_shp_sf <- st_read(dsn = maps_wd, 
                     layer = "sn_100000",
                     stringsAsFactors = FALSE)
## Reading layer `sn_100000' from data source `E:\PopDynIZW Dropbox\Steph Kramer\_GitHub\_DataBorneo\BaseMaps' using driver `ESRI Shapefile'
## Simple feature collection with 396 features and 4 fields
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: 108.9454 ymin: -3.787917 xmax: 118.7879 ymax: 6.464583
## geographic CRS: WGS 84

Load point df and convert to spatial object

These are our observations of the species.

### The spatial point data (species records) (repetition) ###

# filename
spec_pt_filename <- paste(recs_wd,'/','MyNewSpecies.csv', sep='')
spec_pt_filename
## [1] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/_DataBorneo/BaseRecords/MyNewSpecies.csv"
# you can play also with other files
#spec_pt_filename <- paste(recs_wd,'/','DHOsim.csv',sep='')
#spec_pt_filename <- paste(recs_wd,'/','PPLsim.csv',sep='')

# read the file
sp_recs <- read.csv(file = spec_pt_filename, header=TRUE, sep=',')

#convert it to spatial object (sf here)
sp_recs_sf <- st_as_sf(x = sp_recs, 
                       coords = c("long","lat"), # columns  for the coordinates
                       crs = 4326, # define crs, 4326 is the EPSG code
                       sf_column_name = "geometry",
                       remove=F) # sf needs a geometry column and you have to name it

# load a second species 
river_pt_filename <- paste(recs_wd,'/','RIVERsim.csv',sep='')
river_recs        <- read.csv(file = river_pt_filename, header=TRUE, sep=',')
river_recs_sf     <- st_as_sf(x = river_recs,
                          coords = c("long", "lat"), 
                          crs = 4326, 
                          sf_column_name = "geometry")

Plot

### Plot to get an impression ###

plot(ras42, col=grey.colors(20))
plot(PA_shp,border='green', lwd=1.8, add=T)
plot(sp_recs_sf$geometry, pch= '*',cex=1,col='deeppink',add=T)
text(112, 5, 'Starting with SDMs', cex=2, col= 'red')
plot(River_shp[,3], col='dodgerblue4', add=T)

Preconditions

Checking for multicollinearity

# workaround for slow computers. First, aggregate the 1 km² resolution into 50*50 km  
agg_pred <- aggregate(x=predictors,fact=50,FUN=mean)
plot(agg_pred)

Pairs plot

raster::pairs(agg_pred, method = 'spearman')

# Plot the differences (residuals) between the rasters:
diff <- raster::corLocal(x = agg_pred[[1]], y = agg_pred[[2]], method = 'spearman')
plot(diff)

# not_run on slow computers:
#pairs(predictors)

NB: I am a big fan of the pairs()-plot because you actually see the distribution of the data, but you can do really nice plots with the {corrplot}-package:
https://cran.r-project.org/web/packages/corrplot/vignettes/corrplot-intro.html

Sampling bias correction

We are now creating a very rudimentary bias correction file. Let’s assume that doing field work in the province of Sabah, Malaysia, is easy and the field sampling has concentrated on these areas. further, it is easier to sample in lowlands and flat slopes. Usually, protected areas (PAs) are also well studied areas, where the sampling intensity is high. And let’s assume that it was not possible to get a research permit for Indonesia and that therefore very few field studies had been conducted here. Please note: this is a completely fictive example! You should have the information about the sampling bias.

That is, we create a file that is representing these differences in sampling effort. this is easiest to do in raster manipulation. We first create a raster from our shapefile with the country outlines.

head(Borneo_shp_sf)
## Simple feature collection with 6 features and 7 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 115.8098 ymin: 4.27081 xmax: 119.2697 ymax: 7.013332
## geographic CRS: WGS 84
##   ID_0 ISO   NAME_0 NAME_1     NAME_2 Shape_Leng Shape_Area
## 1  158 MYS Malaysia  Sabah Lahad Datu  6.7475367 0.62106012
## 2  158 MYS Malaysia  Sabah      Papar  1.6727672 0.10065591
## 3  158 MYS Malaysia  Sabah  Penampang  0.9548857 0.03951545
## 4  158 MYS Malaysia  Sabah Pensiangan  4.0257859 0.49725547
## 5  158 MYS Malaysia  Sabah      Pitas  2.9223456 0.11955352
## 6  158 MYS Malaysia  Sabah      Ranau  2.4815704 0.24027591
##                         geometry
## 1 MULTIPOLYGON (((118.2299 4....
## 2 MULTIPOLYGON (((116.0144 5....
## 3 MULTIPOLYGON (((116.0477 5....
## 4 MULTIPOLYGON (((116.5357 5....
## 5 MULTIPOLYGON (((117.2762 6....
## 6 MULTIPOLYGON (((116.9182 6....
plot(Borneo_shp_sf[,3]) # column 3 is NAME_0

We create a raster that gets the value of 1 for Malaysia (and Brunei) and a value of 0.01 (MaxEnt cannot deal with 0 -> see lecture) for the rest. We use one of our rasters as mask for cell size resolution and cropping extent:

### *Rasterize* the shapefiles
## command rasterize() 
## field = 1 means all raster cells get value = 1, other cell values are set to 0
## background sets all other cells to a chosen number


Mal_ras <- rasterize(x=Borneo_shp_sf[Borneo_shp_sf$NAME_0 %in%
                       c("Brunei","Malaysia"),], y=ras24, field=1,
                       background=0.1)

plot(Mal_ras)

# almost same as:
# NotInd_ras <- rasterize(x=Borneo_shp_sf[Borneo_shp_sf$NAME_0 !=
#                           "Indonesia",], y=ras24, field=1, background=0.01)

We do the same only for the province of Sabah and the Protected Areas

Sab_ras <- rasterize(x=Borneo_shp_sf[Borneo_shp_sf$NAME_1 == 'Sabah',],
                     y=ras24, field=1, background=0.1)

PA_ras <- rasterize(x=PA_shp, y=ras24, field=1, background=0.1)


### Crosscheck with plot
par(mfrow=c(1,2))
plot(Sab_ras)
plot(PA_ras) #n.b.: only 1 and small values, no NA

par(mfrow=c(1,1))

And we do the same for all elevations below 300 m:

ras24_300 <- ras24 <= 300
plot(ras24_300)

Now we add them all up. This means, we rank the sampling intensity, i.e. sampling in protected areas in Sabah below 300 m gets the highest value

bias_tmp <- ras24_300 + Sab_ras + Mal_ras + PA_ras
bias_tmp
## class      : RasterLayer 
## dimensions : 1425, 1423, 2027775  (nrow, ncol, ncell)
## resolution : 0.008333334, 0.008333334  (x, y)
## extent     : 108.3341, 120.1925, -4.374013, 7.500988  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 0.3, 4  (min, max)
plot(bias_tmp)

Let’s standardize between 0 and 1 (and then set the zeros to 0.01 again)

# get the maximum value in the layer, standardize it and round to two decimals
maxval <- max(values(bias_tmp),na.rm=T)
bias_tmp2 <- bias_tmp/maxval
bias_tmp3 <- round(bias_tmp2, digits=2)

## same as in one go:   
bias1 <- round(bias_tmp/max(values(bias_tmp),na.rm=T),digits=2)

table(values(bias1)) 
## 
##   0.08    0.3   0.33   0.52   0.55   0.75   0.78      1 
##  89415 107956 423878  39157 146727   5611  49196   4634
# table(bias1@data@values) # is doing the same

### in case of having 0 somewhere:
bias1[values(bias1 == 0)] <- 0.01 # because MaxEnt
                                  # does not take 0!
table(values(bias1))
## 
##   0.08    0.3   0.33   0.52   0.55   0.75   0.78      1 
##  89415 107956 423878  39157 146727   5611  49196   4634

Now compare the rasters - do they have the same extent, cell size etc?

compareRaster(ras24,bias1) # the same?
## [1] TRUE

check for the length of the NODATA values:

# Are they really the same???
length(which(!is.na(values(bias1))))
## [1] 866574
length(which(!is.na(values(ras24))))
## [1] 866574

Nice! Now plot it:

### Plot the bias file
plot(bias1, col = viridis(7)) #or col = grey.colors(7)

Save the bias file

## save into Output folder, which is your working directory: 
## check with getwd() if you're not sure of wd

## in case you want it aggregated
# bias_agg <- aggregate(bias1, fact=50,FUN=mean)
# writeRaster(bias_agg, filename="bias_agg.asc", datatype='ascii',
#            overwrite=TRUE,NAflag=-9999)

## since I used setwd() to the output folder, it is automatically stored there,
## if not, paste the path to the filename, i.e. 
## filename=paste0(output_wd,"/bias_2022.asc")
writeRaster(bias1, filename="bias_2022.asc", datatype='ascii',
            overwrite=TRUE,NAflag=-9999)

Now, create folders for the MaxEnt output by hand or do it with the following code. We want a subfolder in our output-folder, called MaxEntRes, with two subfolders called ‘bias’ and ‘nobias’.

#dir.create(path = file.path(output_wd, "MaxEntRes"), showWarnings = FALSE)
#dir.create(path = file.path(output_wd, "/MaxEntRes", "nobias"),showWarnings = FALSE)
#dir.create(path = file.path(output_wd, "/MaxEntRes", "withbias"),showWarnings = FALSE)

Before we work with MaxEnt, please inspect the header of your bias file by opening it in a normal text editor. Mine looks like:

NCOLS 1423
NROWS 1425
XLLCORNER 108.334122887
YLLCORNER -4.374012999
CELLSIZE 0.0083333337997189
NODATA_value -9999

If this throws an error in MaxEnt, please change the
CELLSIZE from 0.0083333337997189 to 0.0083333338 by hand! That is, the header should read:

ncols 1423
nrows 1425
xllcorner 108.33412288693
yllcorner -4.3740129986359
cellsize 0.0083333338
NODATA_value -9999

(Please, without ‘
’ if you use copy-paste from my .rmd course script…..)

#### meanwhile clicking Maxent.... 
# (or: do it with R package dismo, but this is another story
#  which will be taught another time....)

Accessing MaxEnt results

List the results

# Have a look at the MaxEnt layers produced
# MaxEnt saves the mean of x repetitions in _avg.asc

# The script now only works with the folder structure
# MaxEntRes/nobias
# MaxEntres/withbias

# new argument 'recursive' means, that also subfolders are checked!
infiles <- list.files(path=paste(output_wd,'/MaxEntRes',
                                 sep=''),pattern='_avg.asc$',
                      full.names=TRUE,recursive=TRUE )
infiles
## [1] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/output/MaxEntRes/nobias/river_dummy_avg.asc"             
## [2] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/output/MaxEntRes/nobias/river_dummy_FutureMaps_avg.asc"  
## [3] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/output/MaxEntRes/withbias/river_dummy_avg.asc"           
## [4] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/output/MaxEntRes/withbias/river_dummy_FutureMaps_avg.asc"
# example
#[1] "c:/Users/kramer/Dropbox/.../MaxEntRes
#     /nobias/river_dummy_avg.asc"
#[2] "c:/Users/kramer/Dropbox/.../MaxEntRes
#     /nobias/river_dummy_FutureMaps_avg.asc"
#[3] "c:/Users/kramer/Dropbox/.../MaxEntRes
#     /withbias/river_dummy_avg.asc"
#[4] "c:/Users/kramer/Dropbox/.../MaxEntRes
#     /withbias/river_dummy_FutureMaps_avg.asc"

Stack results file and make usual plot

## IF you have a slow computer, 
## please use the workaround with function aggregate() from above:
# me_stack <- aggregate(stack(infiles[c(1:4)]),fact= 50, FUN = mean )

me_stack <- stack(infiles[c(1:4)])

# name sequence according to infiles list above
names(me_stack) <- c('curr_noBias','fut_noBias',
                     'curr_withBias','fut_withBias')

#par(mar=c(4,4,3,3), oma=c(2,2,2,2))

plot(me_stack,col=viridis(100)) # note: I might use a different example than you did in MaxEnt

Check the range of the values

boxplot(me_stack,layers = c(1,3,2,4),notch=T,outline=F)

How would you interpret the results?

Threshold selection

You should have results for the runs with a bias file and without:

# Use maxentResults.csv to access threshold value
# column: 10th percentile training presence logistic threshold

me_res <- list.files(path=paste(output_wd,'/MaxEntRes',sep=''),
                     pattern='maxentResults.csv',
                     full.names=TRUE,recursive=TRUE )
me_res
## [1] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/output/MaxEntRes/nobias/maxentResults.csv"  
## [2] "E:/PopDynIZW Dropbox/Steph Kramer/_GitHub/output/MaxEntRes/withbias/maxentResults.csv"
# example
#[1] "c:/Users/kramer/Dropbox/.../MaxEntRes
#     /nobias/maxentResults.csv"
#[2] "c:/Users/kramer/Dropbox/.../MaxEntRes
#     /withbias/maxentResults.csv"

Now we do some elegant trick: we store both files as two objects in a list! Similar to tacking rasters, a bit….

store_res <- lapply(me_res,FUN=read.csv) #store as list obj.
## Check the following out:
# head(store_res) # list object with 2 slots
# str(store_res)
# store_res[[1]]

We are searching for the column with the threshold.

## Extract the values at column '10th percentile' at position nrep +1
## (= nrow command) to access avg value of all repetitions at the lowest row

zename<-'X10.percentile.training.presence.logistic.threshold' # check if the column name is the same !!!!!
# Achtung! 
# In MaxEnt 3.4.1 heisst Spalte 'X10.percentile.training.presence.Logistic.threshold'
zename
## [1] "X10.percentile.training.presence.logistic.threshold"

Now check which column number that is….

zecol <-which(colnames(store_res[[1]]) == zename) 
zecol
## [1] 73

… and extract the threshold value :

## a little head twister:
## access each element of your list, which is a whole dataset with store_res[[1]]
## and inside this list object = data.frame, access the element of
## the last row (nrow), and the column zecol] containing the threshold value
t_noBias   <- store_res[[1]][nrow(store_res[[1]]),zecol]
t_withBias <- store_res[[2]][nrow(store_res[[2]]),zecol]

t_noBias #the thresholds
## [1] 0.1904
t_withBias
## [1] 0.3014

Now apply the thresholds to your MaxEnt output maps of relative probability. For that, store the thresholds in a vector, and double them for bias/ nobias, because the threshold stays the same:

all_bias <- rep(c(t_noBias,t_withBias), each=2)
all_bias
## [1] 0.1904 0.1904 0.3014 0.3014

Now set the thresholds to the stack. If you don’t remember the stack (me_stack), plot it again.Test different thresholds and see how the predictions change:

binary_thresh <- me_stack >= all_bias
plot(binary_thresh)

test05 <- me_stack >= 0.5 # check result with fixed threshold of default cut-off
plot(test05)

#testSeq <- me_stack >= seq(0.1, 0.7, length = 4)
#plot(testSeq)

Analyses

Now let’s check how global change affects the protected areas. Does the suitability of these areas increase or decrease in the future?

# Extract relative probability values inside protected areas PAs
# check if overlay correct - only works when you have NOT aggregated the maps
compareRaster(PA_ras,me_stack[[1]]) 
## [1] TRUE
ex <- which(values(PA_ras)==1)
# ex # gives Pointer to elements/ Index
ex_stack <- extract(x=me_stack, y=ex)
head(ex_stack)
##      curr_noBias fut_noBias curr_withBias fut_withBias
## [1,]          NA         NA            NA           NA
## [2,]   0.1375010   0.185058     0.1175400    0.0772955
## [3,]   0.0776266   0.127964     0.0860960    0.0611624
## [4,]   0.2222820   0.117761     0.1645070    0.0576882
## [5,]   0.0677608   0.106137     0.0770865    0.0531054
## [6,]   0.5384140   0.337645     0.4211340    0.1714620

Plot the difference of habitat suitability in scenarios and PAs

boxplot(ex_stack,na.rm=T)

Calculate no. of cells (=area) of suitable habitat in PAs

ex_binary <- extract(x=binary_thresh, y=ex)
head(ex_binary)
##      curr_noBias fut_noBias curr_withBias fut_withBias
## [1,]          NA         NA            NA           NA
## [2,]       FALSE      FALSE         FALSE        FALSE
## [3,]       FALSE      FALSE         FALSE        FALSE
## [4,]        TRUE      FALSE         FALSE        FALSE
## [5,]       FALSE      FALSE         FALSE        FALSE
## [6,]        TRUE       TRUE          TRUE        FALSE
## In the following, we are summing the cells with 'TRUE'
ex_binary[1:2,]; colSums(ex_binary,na.rm=T) # 2 cmnds in line with ; !
##      curr_noBias fut_noBias curr_withBias fut_withBias
## [1,]          NA         NA            NA           NA
## [2,]       FALSE      FALSE         FALSE        FALSE
##   curr_noBias    fut_noBias curr_withBias  fut_withBias 
##         43285         63905         43430         39989

Interestingly: when correcting for sampling bias, the future scenarios have the lowest number of suitable habitat inside PAs:

curr_noBias fut_noBias curr_withBias fut_withBias
43285 63905 43430 39989

Least cost path for connectivity

#install.packages('gdistance')
library(gdistance)
#?gdistance

# Set start and end points: centers of the first two PA polygons
start <- gCentroid(PA_shp[1,1], byid=FALSE, id = NULL)
end   <- gCentroid(PA_shp[2,1], byid=FALSE, id = NULL)

# Several commands in one line when separated by ';':
plot(me_stack[[1]]); points(start); points(end)

# Clip raster to save computation time
cr_extent <- c(114,117,3,6.5)
me_cr <- crop(x=me_stack[[1]],y=cr_extent)
plot(me_cr); points(start,pch=15); points(end,pch=15)

### Necessary calculations
# Calculate the transition layer
trans <- transition(x = me_cr,
                    transitionFunction = mean,
                    directions = 4,
                    symm = F)
class(trans)
## [1] "TransitionLayer"
## attr(,"package")
## [1] "gdistance"
# Calculate the shortest weighted connection
sPath <- shortestPath(trans, start, end,
                      output="SpatialLines")

crs(sPath) <- "+proj=longlat +datum=WGS84 +no_defs"


# Calculate the length of the path
costDistance(trans, start, end) #units?
##         1
## 1 1005.26
# Make a plot
plot(me_cr,col=viridis(100)); points(start,pch=15); points(end,pch=15)
lines(sPath,lwd=1, col= 'red'); plot(PA_shp, border='white', lwd=0.5,add=T)